Hierarchical model for the scale-dependent velocity 

of seismic waves 



J. TWORZYDLO^'^ ^ C. W. J. BEENAKKER^ 

""Instituut-Lorentz, Universiteit Leiden, P.O. Box 9506, 2300 RA Leiden, The Netherlands 
'^Institute of Theoretical Physics, Warsaw University, Hoza 69, 00-681 Warszawa, Poland 

Elastic waves of short wavelength propagating through the upper layer of the 
Earth appear to move faster at large separations of source and receiver than 
at short separations. This scale dependent velocity is a manifestation of Fer- 
mat's principle of least time in a medium with random velocity fluctuations. 
Existing perturbation theories predict a linear increase of the velocity shift with 
increasing separation, and cannot describe the saturation of the velocity shift 
at large separations that is seen in computer simulations. Here we show that 
this long-standing problem in seismology can be solved using a model developed 
originally in the context of polymer physics. We find that the saturation veloc- 
ity scales with the four-third power of the root-mean-square amplitude of the 
velocity fluctuations, in good agreement with the computer simulations. 

Seismologists probe the internal structure of the Earth by recording the arrival times of 
waves created by a distant earthquake or explosionB. Systematic differences between studies 
based on long and short wavelengths A have been explained! in terms of a scale dependence 
of the velocity at short wavelengths. The velocity obtained by dividing the separation L of 
source and receiver by the travel time T increases with increasing L, because — following 
Fermat's principle — the wave seeks out the fastest path through the medium (see Fig. p. 
This search for an optimal path is more effective for large separations, hence the apparent 
increase in velocity on long length scales. It is a short-wavelength effect, since Fermat's 
principle breaks down if the width ^/LX of the first Fresnel zone becomes comparable to the 
size a of the heterogeneities. The scale-dependent velocity of seismic waves was noted by 
Wielandt more than a decade agoi, and has been studied extensively by geophysicist^l^S. 

A rather complete solution of the problem for small L was given by Roth, Miiller, and 
Sniederi, by means of a perturbation expansion around the straight path. The velocity shift 
Sv = vq{1 — vqT/L) (with vo the velocity along the straight path) was averaged over spatially 
fluctuating velocity perturbations with a Gaussian correlation function (having correlation 
length a and variance s^Vq, with e <C 1). It was found that {5v) ~ voS^L/a increases linearly 
with L. Clearly, this increase in velocity can not continue indefinitely. The perturbation 
theory should break down when the root-mean-square deviation 5x ~ eaiLjcifl'^ of the 
fastest path from the straight path becomes comparable to a. Numerical simulationsi'i0 
show that the velocity shift saturates on length scales greater than the critical length L^. ~ 
^Qj, validity of perturbation theory. A theory for this saturation does not exist. 
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FIG. 1. Illustration of the scale dependent velocity. Two rays are shown of short wavelength 
waves emitted from a source S and recorded at two receivers R, R' . The shaded areas indicate 
regions of slow propagation. Each ray follows the path of least time from source to receiver, in 
accordance with Fermat's principle. The longer trajectory seeks out the fastest path more efficiently 
than the shorter one, hence the apparent increase in velocity. Perturbation theory breaks down 
when the deviation of the ray from the straight path becomes comparable with the characteristic 
size of the heterogeneities. 

It is the purpose of this article to present a non-perturbative theory for this seismological 
problem, by making the analogy with a problem from polymer physics. The problem of 
the velocity shift in a random medium belongs to the class of optimal path problems that 
has a formal equivalence to the directed polymer problemlll. The mapping between these 
two problems relates a wave propagating through a medium with velocity fluctuations to 
a polymer moving in a medium with fluctuations in pinning energy. The travel time of 
the wave between source and receiver corresponds to the energy of the polymer with fixed 
end points. At zero temperature the configuration of the polymer corresponds to the path 
selected by Fermat's principle. (The restriction to directed polymers, those which do not 
turn backwards, becomes important for higher temperatures.) There exists a simple solvable 
model for directed polymers, due to Derrida and Griffiths^, that has remained unnoticed in 
the seismological context. Using that model we can go beyond the breakdown of perturbation 
theory and describe the saturation of the velocity shift on large length scales. 

Hierarchical model 

We follow a recursive procedure, by which the probability distribution of travel times is 
constructed at larger and larger distances, starting from the perturbative result at short 
distances. At each iteration we compare travel times from source to receiver along two 
branches, choosing the smallest time. A branch consists of two bonds, each bond representing 
the length scale of the previous step. This recursive procedure produces the lattice of Fig. 
0, called a hierarchical latticeS. The lattice in this example represents a two-dimensional 
system, since at each step the length is doubled while the number of bonds is increased 
by a factor of four. For the three-dimensional version one would compare four branches at 
each step (each branch containing two bonds), so that the total number of bonds would 
grow as the third power of the length. Since most of the simulations have been done for 
two-dimensional systems, we will consider that case in what follows. 

To cast this procedure in the form of a recursion relation, we denote by Pk{T) the distribution 
of travel times at step k. One branch, consisting of two bonds in series, has travel time 
distribution 
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FIG. 2. First three steps of the recursive construction of the hierarchical lattice. 



qk{T)= / dT'pk{T')pk{T-T'), (1) 
Jo 

assuming that different bonds have uncorrelated distributions. To get the probabihty dis- 
tribution at step A; + 1 we compare travel times of two branches, 

/■oo /"OO 

Pk+i{T)= dT' dT"5{T-mm{T'X))qk{r)qk{T") 
Jo Jo 

= 2qu{T) / dT'quiT'). (2) 
Jt 

We start the recursion relation at step with the distribution PoiT) calculated from pertur- 
bation theory at length L^. Iteration of eq. (|]) then produces the travel time distribution 
Pk{T) at length L = 2^L^. 

Equation (j^) is a rather complicated non-linear integral equation. Fortunately, it has several 
simplifying propertiesSEl. One can separate out the mean (T)q and standard deviation 
ctq 7^ of the starting probability distribution, by means of the fc-dependent rescaling 
r = {T — 2^ (T)q)/(Jo. The recursion relation @ is invariant under this rescaling, which 
means that we can restrict ourselves to starting distributions Po{t) = <JoPo{aoT + 
having zero mean and unit variance. This is the first simplification. After k iterations the 
mean fhk and standard deviation dk of the rescaled distribution Pkir) yield the mean (T)^ 
and standard deviation ak of the original Pk{T) by means of 

{T)i^ = aoiJik + 2^ (T)o , ak = (Jodk. (3) 



The second simplification is that for large k, the recursion relation for pfc(r) reduces totJ'lla 

Pk+iir) = \apk (|ar + (3dk + (1 - , (4) 

with universal constants a = 1.627 and f3 = 0.647. Under the mapping (^), the mean and 
standard deviation evolve according to 

rhk+i = 2mk - 2(3ak/a, dk+i = 2dk/a. (5) 
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The solution of this simplified recursion relation is 
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(6) 



The coefficients A and B are non-universal, depending on the shape of the starting distri- 
bution po- For a Gaussian po we find A = 0.90, B = 0.95, close to the values A = 1, B = 1 
that would apply if eq. (||) holds down to = 0. For a highly distorted bimodal po we find 
A = 0.71, B = 0.88. We conclude that A and B depend only weakly on the shape of the 
starting distribution. 

Scaling laws 

Given the result (|^) we return to the mean and standard deviation of PkiT) using eq. (^. 
Substituting k = \og2{L/Lc) one finds the large-L scaling laws 
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The mean travel time (T) and standard deviation a scale with L with an exponent p = 
loggO; = 0.702. This scaling exponent has been studied intensively for the directed polymer 
problem^!. 

For the seismic problem the primary interest is not the scaling with L, but the scaling with 
the strength e of the fiuctuations. Perturbation theory! gives the e-dependence at length 



1-Vo (T)o /Lc ~ e Lc/a, VoOq/Lc ~ eJa/Lc, (9) 



where ~ indicates that coefficients of order unity have been omitted. (We will fill these in 
later.) Since Lc — ae~'^/^ (as mentioned in the introduction), we find upon substitution into 
eq. (1^ the scaling of the mean velocity shift at length L ^ Lc. 

{6v) /vo = l- vo (T) /L ~ E^/' [1 + OiLc/LY] . (10) 



The mean velocity shift saturates at a value of order vqe'^^^. The exponent | was anticipated 
in ref. O and is close to the value 1.33 ± 0.01 resulting from simulationJl. 



Comparison with simulations 

For a more quantitative description we need to know the coefficients omitted in eq. (^. 
These are model dependentili'0. To make contact with the simulation^li we consider the 
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FIG. 3. Scale dependence of the velocity, showing the saturation of the mean velocity shift {5v) 
at large separations L of source and receiver. The dashed curve is the result (11) of perturbation 



theory, the solid curves are the non-perturbative results from eq. (13). Data points are results of 
computer simulations at various strengths e of the velocity perturbation (open markers from ref. 
^, filled markers from ref. P). 

case of an incoming plane wave instead of a point source. The perturbation theory for the 
mean velocity shift at length Lc gives 



The variance at length Lc is 




6v^)^ = vy^^V^ 1 1 - V^e'^ ] . (12) 



We quantify the criterion for the breakdown of perturbation theory by Lc = Kae^'^^'^, with 
K = 0.765 our single fit parameter. For the non-universal constants A and B we can use in 
good approximation A = 1, B = 1. The mean velocity shift in the non-perturbative regime 
(L > Lc) is then expressed as: 



+ {Sv)o ■ (13) 




For L < Lc we use the perturbative result ([TT]) (with Lc replaced by L). As shown in Fig. 
3, the agreement with the computer simulations is quite satisfactory, in particular in view 
of the fact that there is a single fit parameter k for all curves. 

Conclusion 

We have presented a non-perturbative theory of the scale-dependent seismic velocity in 
heterogeneous media. The saturation of the velocity shift at large length scales, observed in 
computer simulations, is well described by the hierarchical model — including the e^^^ scaling 
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of the saturation velocity. We have concentrated on the case of two-dimensional propagation 
(for comparison with the simulations), but the e^l"^ scaling holds in three dimensions as well. 
(The coefficients a = 1.74, /3 = 1.30 are different in 3D.) 

Our solution of the seismic problem relies on the mapping onto the problem of directed 
polymers. This mapping holds in the short- wavelength limit, A ^ a? /L. To observe the 
saturation at L ~ thus requires A ^ ae^/'^. For L > a?/\ the velocity shift will decrease 
because the velocitv fluctuations are averaged out over a Fresnel zone. There exists a 
perturbation theoryO for the velocity shift that includes the effects of a finite wavelength. 
It is a challenging problem to see if these effects can be included into the non-perturbative 
hierarchical model as well. 
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